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We analyze differences in dynamics and in properties of the sampled potential energy landscape 
between different equilibrium trajectories, for a system of rigid water molecules interacting with a 
two body potential. On entering in the supercooled region, differences between different realizations 
enhance and survive even when particles have diffused several time their average distance. We 
observe a strong correlation between the mean square displacement of the individual trajectories 
and the average energy of the sampled landscape. 

PACS numbers: 



I. INTRODUCTION 

The Potential Energy Landscape (PEL) 0, 0] is the 
surface generated by the potential energy of a system. 
In the case of system composed by N rigid molecules, 
it is an highly complex surface defined in a 6iV di- 
mensional space. In recent years, numerical studies 
SSHElzLllUQiJInl, boosted by the increased 
numerical resources, have attempted to quantify the sta- 
tistical properties of this surface (for example quantify- 
ing the distribution in energy of specific points of the 
surface, like local minima and saddles) in the attempt to 
develop a thermodynamic description fully based on PEL 
properties. This line of thinking, pioneered by Frank 
Stillinger and his coworkers has been very fruit- 

ful in the study of supercooled liquids, both in equi- 
librium [3 , la, |L3| and in out-of equilibrium condi- 
tions \l4. Il5l [TbL Il7l|. Stillinger's formalism builds on 
the idea that the PEL surface can be partitioned into 
disjoint basins. A basin is unambiguously defined as 
the set of points in configuration space connected to the 
same local minimum — named inherent structure (IS) 
- via a steepest-descendant trajectory. In this respect, 
the PEL's statistical properties entering in the evalua- 
tion of the partition function are the number, shape and 
depth of the PEL basins. For its conceptual simplicity 
and its strict connection with numerical implementation, 
the PEL formalism has become one of the modern tools 
to interpret and analyze simulation data. 

Water, due to its intrinsic interest as liquid of life, has 
been extensively studied in computer simulations. Mod- 
els have been developed which are able to reproduce qual- 
itatively the thermodynamic and dynamic anomalies of 
this liquid B H|. Indeed, water is characterized 
by a line of isobaric density maxima and compressib ility 
minima. Its specific heat C p increases on cooling [2lll22| . 
Recent studies have attempted to connect the thermody- 
namic anomalies to the presence of a line of first order 
transition between two liquid structures, eventually end- 
ing into a second order critical point [23, [24J . Dynamics 
in water also shows anomalies. The diffusion coefficient 
decreases as pressure is increased up to a maximum value 
(approximately 400 MPa[25j ) at room temperature. For 



greater pressures, water dynamics becomes progressively 
slower, as expected for simple liquids. 

Water simulations have been among the first to be 
scrutinized using the methods developed by Stillinger 4, 
II Sill El El EH in the attempt to clarify differ- 
ences in local structures and relations between structural 
properties and dynamics. More recently, extensive stud- 
ies of the PEL properties have been published, provid- 
ing a detailed description of the landsca pe properties 
of several models of water 0, 0, El l33l l34j. Within 
the PEL formalism, attempts have also been presented 
in the direction of connecting thermodynamic proper- 
ties (like the number of explored minima or the num- 
ber of diffusive directions in the PEL) to dynamics 

11 11, 111 in m m m m m m m. The outcome 

of these works, and related studies on different models 
of glass forming liquids, suggests a strong connection 
between dynamics and landscape properties. Still, the 
interpretation schemes require fitting parameters whose 
physical interpretation in often unclear. As an exam- 
ple, it has been shown that ln(D) is consistent with 
a l/(TS con/ ) dependence 0, 1 H M (being S conf 
the configurational entropy) but no understanding of the 
proportionality coefficient in term of PEL properties has 
been reported |45| . 

In this article we report analysis of distinct equilibrium 
trajectories of a system of 216 water molecules interacting 
via a two body potential. The novel aspect is the pos- 
sibility of comparing a large number of independent tra- 
jectories (more than 100 for each state point) for a large 
number of state points. Each of these (previously equili- 
brated) independent trajectories last more than 20 ns, at 
which time — even at the slowest studied temperature — 
the average displacement of the molecule is longer than 
three nearest neighbors. The analysis of this set of data 
shows that, while individually each trajectory appear to 
be to a very good approximation a diffusive trajectory 
(i.e. with a mean square displacement which increases 
linearly with time) different realizations are character- 
ized by apparent values of the diffusion coefficient which 
differ by more than one order of magnitude, a clear ev- 
idence of dynamic heterogeneities. We also show that 
differences in dynamics are strongly coupled with the lo- 
cation of the system on the PEL, providing evidence that 
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dynamic heterogeneities are related to inhomogeneities 
in the local basin energies. For completeness, we briefly 
discuss the distribution of IS energies for this model and 
the temperature and density (p) dependence of the dif- 
fusion coefficient, presenting new data with significantly 
improved statistics. 



II. SIMULATION DETAILS 

We analyze molecular dynamics trajectories of a sys- 
tem composed of 216 rigid water molecules in the (NVE) 
ensemble. Molecules interact via the widely studied two- 
body simple point charge extended (SPC/E) model [l8|| . 
The integration step is 1 ft second, and long range in- 
teractions are taken in account using the reaction field 
method. Dynamics and thermodynamics properties for 
the SPC/E model have been extensively studied in the 
past ^(| in a wide range of temperatures and densi- 
ties. Analysis of the statistical properties of the po- 
tential energy landscape have also been characterized 
0, IE HE HE E3 • Here we analyze a large set of 
data, based on 5 different temperatures in supercooled 
states and nine different densities, from 0.9 to 1.4 g/cm 3 . 
At each state point, we analyze equilibrated trajectories 
of more than 100 independent realizations, to both re- 
duce significantly the numerical error and to estimate 
the self-averaging properties of this model for different T 
and p. Each trajectory covers a time interval of about 
20 ns. Potential energy landscape properties have been 
based on the analysis of the inherent structures, which 
has been calculated using a standard conjugate gradient 
minimization algorithm with I0 -15 tolerance |47j . 



III. IS ENERGIES: THE RANDOM ENERGY 
MODEL 

A quantification of the statistical properties of the po- 
tential energy landscape, i.e. the distribution of basin's 
depth and shapes, for the SPC/E model of water has 
been recently reported j3^|. It has been shown that the 
number Q(ejs)deis of distinct basins of energy depth eis 
between ejg and ejg + deis follows a Gaussian distribu- 
tion, 



tt(eis)dei S = e 



-(e IS -E ) 2 /2<r 2 



V2 



7T(7 



-Cfe/s, (1) 



where e is the total number of distinct basins of the 
PEL for the system of N molecules, Eq is the energy 
scale of the distribution and a 2 is the variance. The co- 
efficients a, E Q and a depend only on the system density. 
The corresponding probability distribution P(e/s,T) of 
sampling an IS of depth eis in equilibrium at tempera- 



ture T is given by 
P(e IS ,T) 



JP(e IS ,T)de IS 



(2) 



where (3 = 1/kT and f V ib(T, eis) is the basin free energy 
[HOI- The hypothesis of a Gaussian form for fi(esj), 
together with the assumption of an eis independence of 
the basin anharmonicities implies: (i) that the T depen- 
dence of the average IS energy < eis(T) > at constant 
volume, is linear as a function of T _1 £IEEEjIE3|, that is: 



<ei S (T) >=A+|, 



(3) 



where the coefficients A and B depend only on p and can 
be expressed in term of landscape properties[5(|, and (ii) 
that the probability distribution P(ejs,T) is Gaussian 
with a T-independent variance jH, |IE HE LlJ HE HE H3] • 
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FIG. 1: Temperature dependence of the average < eis > (T) 
for all investigated densities. 

To provide evidence of this behavior we report in Fig. ^ 
< eis(T) > at all studied densities. In all cases, in the in- 
vestigated T-range, the expected T~ ^dependence holds. 
Nevertheless, two breakdowns of Eq|3]are expected out- 
side the investigated T-region. At high T, due to the 
increasing importance of the anharmonic contribution to 
the basin free energy and at low T due to possible 
presence of non-Gaussian corrections to the eis proba- 
bility distribution (Eq. {Q. The low T region can not be 
numerically studied due to the huge increase of the relax- 
ation times. Fig. [21 shows the p dependence of A and B 
(which have well defined landscape interpretations 50] ). 
As a first approximation (but see Ref. 0] for a more pre- 
cise discussion), A is related to E and B to a 2 . 

The large available data set allows us to evaluate 
P(eis,T) and test the second prediction of the Gaussian 
hypothesis. Fig.^a and^Ja show P(eis, T) respectively 
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FIG. 2: Density dependence of the fitting parameters A and 
B in Eq0 



for p = 0.95 g/cm 3 and p = 1.00 g/cm 3 . Figures show 
that the shape of the distribution is, to a first approxima- 
tion, Gaussian and that the variance of the distribution 
only weakly increases with T, as can be inferred by the 
small change of the heigh of the distribution with T. 

The availability of 100 independent realizations allows 
us to study also the distribution of P(< eis >i,T), de- 
fined as the distribution of the average depth sampled in 
each independent trajectory i in the simulated 20 ns time 
interval. In this case, of course, each distribution is evalu- 
ated only over 100 points. If the length of each trajectory 
is sufficiently long, so that the system is able to sample 
all basins which are statistically relevant at that tem- 
perature, the distribution should be peaked around the 
average of P(eis, T), with a small variance. Fig.^b and 
Fig. 0Jb show that this is the case only at the higher T. 
On cooling, the relaxation time of the system increases 
and, within the time of the simulation, the system retains 
memory of the initial basin. It is important to observe 
that the distribution becomes very asymmetric, develop- 
ing a long tail at low basin energies. This may suggest a 
strong relation between depth of the explored basin and 
dynamics, which can not be attributed to differences in 
thermal energy, since all different realizations have the 
same kinetic energy. In particular, the asymmetry in the 
resulting distribution suggests that configurations start- 
ing from low energy basins do not have time to explore 
phase space sufficiently. It is worth stressing, that similar 
results are observed at all investigated densities. To sup- 
port this hypothesis we show in Fig.[S]the time evolution 
of e/s and of the center of mass mean square displace- 
ment MSDi for two different trajectories at the same T 
and p. More precisely, the mean square displacement of 
the trajectory i is defined as 



MSDi 



1 N 
N l^Vj 



where fy (t) is the location of the center of mass of 
molecule j in the trajectory i at time t. The averaged 
M SD is defined as 

l M t 

MSD = —r V MS* A, (5) 
Alt t— 1 

where Aft is the number of independent trajectories. The 
two chosen trajectories differ by the value of ejs at time 
0. In the studied t range (ns), both configurations sam- 
ple a restricted interval of e/s values. The memory of the 
energy of the starting basins is preserved during the sim- 
ulation time. Comparing Fig. 0-a and Fig. 0-b we note 
that the configuration with large e/s is characterized by 
a larger MSD. 
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FIG. 3: a) Probability distribution of the eis (per molecule) 
at density d — 0.95 g/cm 3 , for different T values, b) Proba- 
bility distribution of < eis > , where each average is calcu- 
lated over one distinct trajectory. Units are arbitrary. Each 
P(eis,T) is calculated using 5000 points, since from each of 
the 100 independent trajectories we evaluated 50 different in- 
herent structures, minimizing at intervals equally spaced in 
time. 



IV. DYNAMIC HETEROGENEITIES 

We next focus the attention on dynamic hetero- 
geneities via a study of the mean square displacement of 
the individual trajectories [3^, |5^- For each state point 
(T, p) we calculate MSDi for each of the Aft trajectories. 
Fig. EJshows MSDi for p = 0.95 g/cm 3 at three selected 
T, but similar results are obtained for all the other stud- 
ied densities. In all cases, no averaging over different 
time origins is performed. In the figure, the time axis 
are chosen in such a way that the average MSD at the 
maximum reported time coincides for all temperatures. 
A comparison of the spreading of the different realiza- 
tions at fixed value of the average MSD, for example at 
the maximum reported time, reveals a clear increase in 
the fluctuations of the different trajectories on lowering 
T. 
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FIG. 4: Same as FigEJfor d = 1.40 g/cm 3 




FIG. 5: eis and the mean square displacement as a function 
of time for two different trajectory at T = 21QK and p = 0.95 
g/cm 3 . 



To better quantify fluctuations in dynamics and the 
role of T, we calculate the variance of the mean square 
displacement of different trajectories. More precisely, we 
evaluate, for each time t 



z2i(MSDi — MSD) 2 



(6) 



Fig. shows the behavior of the variance <J MSD as a 
function of the average MSD, parametrically in t, for 
several different T. This representation, which accounts 
for the intrinsic effect of the slowing down of the dynam- 
ics on cooling by eliminating t, confirms that dynamic 
heterogeneities grow significantly on supercooling. Data 
in Fig.[7|quantify that the spreading of the MSDi values 
is much more enhanced at low T. The T-dependence of 
the slope of the curves shown in Fig 0is shown in Fig. [5] 

It is interesting to compare results reported in Fig. [7| 
and Fig. [SJwith expectations for a Gaussian random walk 
process. For this case, the probability of finding the 




FIG. 6: Individual mean square displacement for 100 inde- 
pendent realizations at three different temperatures. Filled 
symbols indicate the MSD averaged over the different real- 
izations. Note that no average over starting time has been 
performed. 



walker after time t at distance r from the location at 
the time origin is 



P{r 2 ,t)dr 2 



2n<r 2 (t) > 3 / 2 



e 2<rHt)> dr 2 . (7) 



The first and second moment of this distribution are 
given by 



<r z >= r A P{r\t)dr z =<r\t)>, (8) 
'o 
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FIG. 7: Variance o MSD of the MSD as a function of the 
average MSD for p = 0.95 g/cm 3 , for different T. The two 
full lines indicate the two extreme limits provided by Eqs. lilt 
1121 Note that in water the nearest neighbor distance is 0.28 
nm, corresponding to a square displacement of ~ 0.09 nm 2 . 

and 

< r 4 >= / r 4 P{r 2 ,t)dr 2 = £ < r 2 (t) > 2 . (9) 
Jo 3 

The variance 5 of r 2 between different trajectories of 
the walker is 

S 2 =< (r 2 - < r 2 >f >=\< r 2 (t) > 2 . (10) 

Hence, a plot of 5 vs. < r 2 (t) > has, for a single 
Gaussian walker, a slope of yf, an universal value inde- 
pendent on the diffusion constant. 

If the dynamics in simulated water could be repre- 
sented by the dynamics of N = 216 independent walker, 
then omsd should be related to MSD by the relation 

vmsd = ^^=MSD. (11) 

since each MSDi would be the sum of N independent 
gaussian processes (with a reduction of the variance by 
a factor y/~N as compared to the single random walker 
case). Data in Fig. [7| and Fig. [S] shows that this limit is 
approached at high T. 

In an other extreme theoretical case, each realization 
can be considered as a single random gaussian process 
(for example, in the limit of strong correlation between 
all N molecules, or in the limit of one single diffusing 
molecule). In this limit, the expected relation between 
o'msd and MSD would be 

& msd = J\MSD. (12) 

For reference, this limit is also reported in Fig. [7| and 
Fig. E] 



FIG. 8: Slope of a MSD vs. MSD as a function of T for p = 
0.95 g/cm 3 . The two lines indicate the two extreme limits 
provided by Eqs. 1111121 

Data shown in Fig. and Fig. [S] show that a cross-over 
from the behavior of Eq. ^] to the behavior of Eq. IT21 
takes place on supercooling. It will be very interesting 
to study the size, T and t dependence of these effects. 
Moreover, since each of the equilibrium trajectory of the 
216 molecules system can be thought as representative 
part of a large system, the increase of the variance with 
decreasing T provides a strong evidence of growing dy- 
namic heterogeneities in the system. In this respect, this 
set of data, or analogous data for simpler potentials, may 
become a relevant tool for discriminating between differ- 
ent theories of the glass transitions, in particular between 
the ones based on facilitated dynamics ideas [5(| and trap 
models [13, H| . 



V. CORRELATION BETWEEN DIFFUSION 
COEFFICIENT AND e IS 

Bulk dynamic properties of SPC/E water have been 
previously investigated in details |46l l59l l60l l6ll l62L l63j . 
In particular, PEL inspired studies have investigated the 
relation between the T and p dependence of the diffusion 
coefficient D and the number of unstable directions in 
configuration space [H, as well as the relation between 
D and the configurational entropy 4] . Fig. [5] shows the 
average diffusion coefficient D, evaluated from the MSD 
long time limit, for several studied density and temper- 
atures. The present set of data improves the precision 
of previous estimate for the same model |(32|. Despite 
the relatively small T-range investigated in this study, D 
varies over more than three order of magnitudes. It also 
clearly show that, among all the studied densities, D is 
largest around density 1.1 g/cm 3 . Increasing p, dynamics 
slow down due to packing effect, while decreasing p dy- 
namics slow down due to the development of a network 
of hydrogen bonds. 

Most of previous studies on dynamic 

heterogeneities |6^, |6(| have focused on the dy- 
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FIG. 9: Average diffusion coefficient D as a function of T for 
several studied densities. 



namics of different sub-set of molecules in the same 
system, raising sometime the question if the observed 
results were somehow associated with the chosen rules 
for identifying slow and fast particles. It is also fairly 
difficult to correlate properties of the subset of molecules 
to energetic properties, due to the difficulty of separating 
unambiguously the single particle energy contributions. 
In the present approach, heterogeneities are addressed 
globally, as a fluctuation phenomenon. To better clarify 
the connection between dynamic and energetic (static) 
heterogeneities, we correlate the apparent diffusion coef- 
ficient of each trajectory with the corresponding average 
inherent structure. Indeed, it is important to observe 
that, when each trajectory is averaged over different time 
origins, the resulting MSD behaves linearly with t, such 
that an apparent trajectory diffusion coefficient Di can 
be estimated. Differences in Di values between different 
trajectories persist even for time such that molecules 
have diffused over distances larger than two molecular 
diameters (being « 0.3 nm the distance between the 
center of mass of two nearest neighbor molecules). 

Fig. 1 1UI shows a plot of the average D(T) as a function 
< ejs > (T) for one selected isochore. It also show for 
each trajectory i ( at the same T and p), the apparent 
diffusion coefficient Di vs the average sampled < e/s >i- 
While at larger T, self-averaging is well accomplished on 
the time scale of the simulation, at lower T, each trajec- 
tory — even if diffusive over distances of several molecu- 
lar diameters — samples only a small part of the statis- 
tically relevant configurational space, resulting in a large 
spreading in the value of < e/s >i- Interestingly enough, 
such a spreading in < e/s >, is strongly correlated to the 
spreading in the D, values. 

The relation between Di and < e/s >j displayed by the 
data is not very different from the corresponding relation 
D vs < ejs > for the averaged values. A more detailed 
study of Di vs < e/s >i over a smaller grid of temper- 
atures, allowing for overlap of different (£>;— < e/s >;) 
pairs of points could help sorting out the T and the e/s 
roles. Indeed, one could associate each e/s to a precise 
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FIG. 10: Comparison between the average diffusion coefficient 
(line with open circles) and the diffusion coefficient of each 
single trajectory at p — 0.95 g/cm 3 (top panel) and p = 1.40 
g/cm 3 (lower panel). 



value of Z?(e/s) and in analogy with the thermodynamics 
formalism developed by Stillinger, one could attempt to 
calculate D(T) as[ll| 

D(T) = J D(e IS )T(e IS ,T)P(e IS ,T)de I s (13) 

, where ^(e/s,? 1 ) account for the role of T in dynam- 
ics and P(eis,T) is calculated according to Stillinger's 
PEL formalism. If this goal would be reached, the PEL 
formalism would become an exceptionally reach tool not 
only for describing the thermodynamics of supercooled 
liquids but also their dynamics. 



VI. CONCLUSION 

One of the key features in the slowing down of the dy- 
namics of liquids on approaching supercooled states - 
i.e. when dynamics begin to slow down significantly as 
compared to the standard liquid values — is the devel- 
opment of local fluctuations both in static and dynamic 
properties. The development of numerical resources al- 
lows us to start looking carefully into this problem, by 
performing analysis of the fluctuations|54| . These ap- 
proaches, compared to the corresponding studies of the 
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average properties, are complicated by the interplay be- 
tween space and time. The data reported in this article, 
more than providing conclusive answers, hopefully clarify 
the richness of this type of analysis and stimulate further 
studies focusing on the time and space evolution of the 
fluctuations. 

The presented preliminary analysis reported clearly 
show that the role of fluctuations significantly grows 
on cooling. Comparing fluctuations in dynamics at 
fixed mean square displacement (Fig[7|), we have de- 
tected a progressive increase of the fluctuations on cool- 
ing, already in the region of early supercooling, where 
MCT appears to provide a consistent description of the 
dynamics[H3|- Interestingly, self-averaging properties ap- 
pear to set in only after molecules have diffused several 
particle diameters. This finding not only clarify the diffi- 
culty of calculating reliable values for dynamical quanti- 
ties in deep supercooled states but call attention on the 
fact that a complete decorrelation of the system requires, 
at low T, rearrangements which extends much beyond the 
first neighbor shell. 

We have also observed a clear correlation between the 



apparent diffusion coefficient of the individual realiza- 
tions and depth of the sampled PEL (Fig ll(J|) . Again, 
the comparison of different trajectories, all at the same 
temperature, helps in eliminating trivial (but a priory 
unknown) thermal effect, highlighting the connection 
between dynamics and IS energies. In this respect, one 
of the goal of future studies should be the development 
of a (size dependent) dynamical histogram reweighting 
formalism, conceptually similar to the one used to 
calculate the density of states, from which the relation 
D(eis) could be extracted. This would allow us to sort 
out the role of T and e/s in dynamics, and to describe 
in term of PEL properties not only the thermodynamics 
of supercooled liquids but also their dynamics. 
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